/*
 * Reaction.cpp
 *
 *  Created on: 23 May 2011
 *      Author: allan
 */

#include "Reaction.h"

// C++ includes
#include <cmath>
#include <iomanip>
#include <sstream>
#include <tuple>

// Boost includes
#include <boost/foreach.hpp>

// GeoReact includes
#include <Utils/Algorithms.h>
#include <Utils/Assert.h>

ReactionEquation::ReactionEquation()
{
}

ReactionEquation::ReactionEquation(const vector<string>& species, const vector<double>& stoichiometries) :
species(species), stoichiometries(stoichiometries)
{
}

Reaction::Reaction()
{
	SetKinetics(ZeroKinetics());
}

Reaction::Reaction(const ReactionEquation& equation, const Multiphase& multiphase) :
species(equation.species), stoichiometries(equation.stoichiometries),
indexes_species(multiphase[species])
{
	SetKinetics(ZeroKinetics());
}

Reaction::Reaction(const ReactionEquation& equation, const EquilibriumConstant& K, const Multiphase& multiphase) :
species(equation.species), stoichiometries(equation.stoichiometries),
indexes_species(multiphase[species]), K(K)
{
	SetKinetics(ZeroKinetics());
}

void
Reaction::SetEquilibriumConstant(const EquilibriumConstant& K)
{
	this->K = K;
}

const vector<string>&
Reaction::GetSpecies() const
{
	return species;
}

const vector<double>&
Reaction::GetStoichiometries() const
{
	return stoichiometries;
}

const vector<Index>&
Reaction::GetSpeciesIndexes() const
{
	return indexes_species;
}

const EquilibriumConstant&
Reaction::GetEquilibriumConstant() const
{
	return K;
}

const Kinetics&
Reaction::GetKinetics() const
{
	return *kinetics;
}

const double
Reaction::Stoichiometry(const string& s) const
{
	const unsigned idx = IndexOf(s, species);

	return (idx < species.size()) ? stoichiometries[idx] : 0.0;
}

const double
Reaction::EquilibriumConst(double T, double P) const
{
	return K(T, P);
}

const double
Reaction::ReactionQuotient(const VectorXd& a) const
{
	double Qr = 1.0;

	for(unsigned i = 0; i < species.size(); ++i)
	{
		const double ai = a[indexes_species[i]];
		const double vi = stoichiometries[i];

		Qr *= pow(ai, vi);
	}

	return Qr;
}

const double
Reaction::KineticRate(double T, double P, const VectorXd& n, const VectorXd& a) const
{
	return kinetics->KineticRate(*this, T, P, n, a);
}

ostream&
operator<<(ostream& out, const Reaction& reaction)
{
	const vector<double>& stoichiometries = reaction.GetStoichiometries();
	const vector<string>& species = reaction.GetSpecies();

	string lhs, rhs;

	for(unsigned i = 0; i < species.size(); ++i)
	{
		stringstream stream;

		if(stoichiometries[i] > 0)
		{
			if(stoichiometries[i] != 1.0)
				stream << stoichiometries[i];
			stream << species[i] << " + ";
			rhs.append(stream.str());
		}
		else
		{
			if(stoichiometries[i] != -1.0)
				stream << -stoichiometries[i];
			stream << species[i] << " + ";
			lhs.append(stream.str());
		}
	}

	// Erase the leading string "+ "
	lhs.erase(lhs.size() - 3, 3);
	rhs.erase(rhs.size() - 3, 3);

	out << lhs << " = " << rhs;

	return out;
}

